#include "tspline.h"

TSpline::TSpline()
{
     x=a=b=c=d=NULL;
     n=0;
}

TSpline::~TSpline()
{
    delete[] a;
    delete[] b;
    delete[] c;
    delete[] d;
    delete[] x;

    x=a=b=c=d=NULL;
}

double TSpline::getValue(double value)
{
    int i = getInterval(value);
    double razn = value-x[i];
    return a[i]+b[i]*razn+c[i]/2.0*razn*razn + d[i]/6.0 *razn*razn*razn;
}

double TSpline::getDerivativeValue(double value)
{
    int i = getInterval(value);
    double razn = value-x[i];
    return b[i]+c[i]*razn + d[i]/2.0 *razn*razn;
}

double TSpline::getSecondDerivativeValue(double value)
{
    int i = getInterval(value);
    double razn = value-x[i];
    return c[i] + d[i]*razn;
}

void TSpline::createSpline(const QVector<double> &x_, const QVector<double> &f, const int n_, const double S0, const double Sn)
{
    n=n_;
    delete[] a;
    delete[] b;
    delete[] c;
    delete[] d;
    delete[] x;
    a = new double [n+1];
    b = new double [n+1];
    c = new double [n+1];
    d = new double [n+1];
    x = new double [n+1];
    memcpy(x,&x_[0],x_.size()*sizeof(x[0]));

    double* alpha = new double [n+1];
    double* beta = new double [n+1];
    memset(a,0,(n+1)*sizeof(double));
    memset(b,0,(n+1)*sizeof(double));
    memset(c,0,(n+1)*sizeof(double));
    memset(d,0,(n+1)*sizeof(double));
    memset(alpha,0,(n+1)*sizeof(double));
    memset(beta,0,(n+1)*sizeof(double));

    double h =  (x[n]-x[0]) / n;
    double fi = 0.0;

    alpha[1]=0.0;
    beta[1]=S0;
    a[0]=f[0];
    a[n]=f[n];
    c[0]=S0;
    c[n]=Sn;

    for (int i=1; i<n; i++)
    {
        a[i]=f[i];
        alpha[i+1] = -1.0/(4+alpha[i]);
        fi = 6.0*(f[i-1]-2.0*f[i]+f[i+1]) / h;
        beta[i+1] = - (beta[i] - fi/h )/(4+alpha[i]);
    }

    for (int i=n-1; i>=0; i--)
    {
        c[i]=alpha[i+1]*c[i+1]+beta[i+1];
        d[i+1] = (c[i+1]-c[i]) / h;
    }

    for (int i=0; i<n; i++)
    {
        b[i+1] = (a[i+1]-a[i])/h + c[i+1]/2.0*h - (c[i+1]-c[i])/6.0*h;
    }
    delete[] alpha;
    delete[] beta;
}

int TSpline::getInterval(double value)
{
    int i = 1;
    while (value>x[i])
    {
        i++;
        if (i==n)
            break;
    }
    return i;
}
